function rest=checkou(alpha,m,y0,beta,pathnum)


intBMsq=zeros(100,1);
time=linspace(0,1,100)';
time2=linspace(0,1,10000);
for pn=1:pathnum
    disp(pn);
    getpath=ouprocess(alpha,m,y0,beta,10000); % get a ou process, or bm
    BMsq=getpath.^2; %bm's square
    %BMsq=sin(time2).^2; %sin(s)
    %BMsq=exp(time2).^2;
    %BMsq=getpath;
    for jj=1:100
        intBMsq(jj)=intBMsq(jj)+midptrule(BMsq(1:(jj-1)*101+1),1/9999)^2;
    end
end

rest=intBMsq/pathnum;
hold on
%plot(time,1/2*time.^2)     %case 1 E int_0^t Ws^2ds
%plot(time,(1/2*time-1/4*(sin(2*time))).^2); %case 2 E int0t sin^2sds
%plot(time,1/4*(exp(2*time)-1).^2) %case 3 E int0t e^2sds
plot(time,5/8*time.^4)
plot(time,rest,'r-')




% t=1;
% mcmean=0; mcvar=0;
% for ii=1:pathnum
%     currentpath=ouprocess(alpha,m,y0,beta,10000);
%     mcmean=mcmean+currentpath(10000);
%     mcvar =mcvar +currentpath(10000)^2;
% end
% mcmean=mcmean/pathnum; mcvar=mcvar/pathnum-mcmean^2;
% exmean=y0*exp(-alpha*t)+m*(1-exp(-alpha*t));
% if alpha ~=0
% exvar =beta^2/2/alpha*(1-exp(-2*alpha*t));
% else
% exvar =1;
% end
% meanerr=abs(exmean-mcmean);
% varerr =abs(exvar -mcvar );

